# price update formula
price.update <- function(price,nu) {
  #market shares
  pim <- prob.mat(price)
  # apply the first order condition for profit-max. # c -(s+r)/(-d)
  p.star <- mc - (colMeans(pim)+cannib(pim,price))/ownderiv(pim)
  # weight the nu price
return( as.vector(nu*p.star +(1-nu)*price) )
}
#prob.mat creates the whole prob matrix 
prob.mat <- function(p) {
  # matrix of indirect utilties: row= household, col = variety
  ubar <- b.h %*% t(x) - alpha.h %*% t(p) +rep(1,n.hh) %*% t(xi)
  # MATRIX [n.hh X n.models] of household-level probabilities of variety j
  exp(ubar)/(U0+rowSums(exp(ubar)))
}
#an n.models length vector of derivatives of OWN share wrt OWN price
ownderiv<- function(prm) colMeans(-alpha.h*prm*(1-prm))
#an n.models length vector of derivatives of OTHER CO-OWNED shares wrt OWN price
cannib <- function(prm,price)  {
 # D.cross <- build.DeltaX(prm)
  D.cross <- crossDeltaMLOG(MM=co.own,HH=prm,A=alpha.h)
   return(D.cross %*% (price-mc))
}
#matrix of cross price derivatives, changing price of m (row), impact on j(cols)
build.DeltaX <- function(prm) {
  Delta <- matrix(0,nrow=n.models,ncol=n.models)
  for(i in 1:n.models) {
    for(j in 1:n.models) {
      if(co.own[i,j]) Delta[i,j] = mean(alpha.h*prm[,i]*prm[,j])  
    }
  }
  return(Delta)
}
#n.b. we don't use omega() function in the build.DeltaX as we need to do in MCES.
# also note that our Delta has zeros for own price terms (unlike BLP)
#prm is the vector of household level probabilities
ownelas<- function(prm,share,price) {  # note: ownelas <0 , so it is -epsilon_m in the paper
  #average partial effect of price on market share
  ape <- colMeans(-alpha.h*prm*(1-prm))
  # return a n.models length vector of own-price elasticities
  return(ape*price/share)
}
#singleproduct  super elas
superelas<- function(prm,share,price) {
  #ownprice elas
  ape <- colMeans(-alpha.h*prm*(1-prm)) #average partial effect of price on market share
  eps <- ape*price/share # n.models length vector of own-price elasticities
  #ownderiv elas
  num <- colMeans(alpha.h^2 * prm * (1-prm) * (1-2*prm))
  ownderivelas <- price * (num/ape)
  # return a n.models length vector of superelasticities
  return(1-eps+ownderivelas) # 1+eps in paper
}
#
#variety-level average income of consumers   
avg.inc <- function(prm) {
  num <- colMeans(prm*y.h) 
  denom <- colMeans(prm)
  num/denom
}  
#an n.hh length vector of average elasticities
hhelas<- function(p) rowMeans( (-alpha.h %*% t(p)) * (1-prob.mat(p)) )
